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Abstract 

We propose a general approach to the question of how biological rhythms 

spontaneously self-regulate, based on the concept of "stochastic feedback". 

We illustrate this approach by considering the neuroautonomic regulation 

of the heart rate. The model generates complex dynamics and successfully 

accounts for key characteristics of cardiac variability, including the 1// power 

spectrum, the functional form and scaling of the distribution of variations, and 

correlations in the Fourier phases. Our results suggest that in healthy systems 

the control mechanisms operate to drive the system away from extreme values 

while not allowing it to settle down to a constant output. 
PACS numbers: 02.50.Ey, 05.40, 05.45 
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The principle of homeostasis asserts that biological systems seek to maintain a constant 
output after perturbation |]J. Recent evidence, however, indicates that healthy systems even 
at rest display highly irregular dynamics [Q]. Here, we address the paradox of how to rec- 
oncile homeostatic control and complex variability The concept of dynamic equilibrium or 
homeostasis |]J led to the proposal that physiological variables, such as the cardiac interbeat 
interval T(n), maintain an approximately constant value in spite of continual perturbations. 
Thus one can write in general 

r(n) = T + r], (1) 

where r is the "preferred level" for the interbeat interval and 7] is a white noise with strength 
a, defined as the standard deviation of r\. 

We first re-state this problem in the language of random walks. The time evolution of an 
uncorrelated and unbiased random walk is expressed by the equation r(n+l)— t(ji) = 77. The 
deviation from the initial level increases as n 1 / 2 , so an uncorrelated and unbiased random 
walk does not preserve homeostasis (Fig. |I]a). To maintain a constant level, there must be 
a bias in the random walk ||, 

r(n + 1) - r(n) = I(n) , (2a) 

with 

( w (1 + rj) , if r(n) < r , 
I(n) = (2b) 
[-w (I+77) , if r(n) > r . 

The weight w is the strength of the feedback input biasing the walker to return to its 

preferred level To- 

We find that Eq. (0) does not reproduce the statistical properties of the empirical data 
(Fig. |l]b). We therefore generalize Eq. (0) to include several inputs Ik (k — 0, 1, • • • , m), with 
different preferred levels Tfe, which compete in biasing the walker, and Eq. (0a) becomes 

m 

r(n + l)- r(n) = J2h(n), (3a) 

k=0 
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where 




w k (1 + 77) , if r(n) < r k , 



(3b) 



(I + 77) , if r(n) > r fc . 



From a biological point of view, it is clear that the preferred levels r k of the inputs I k 
cannot remain constant in time, for otherwise the system would not be able to respond to 
varying external stimuli. We assume that each preferred interval r k is a random function 
of time, with values correlated over a time scale T^ ck . We next coarse grain the system 
and choose Tfc(n) to be a random step-like function constrained to have values within a 
certain interval and with the length of the steps drawn from a distribution with an average 
value Ti oc k (Fig. |l]c). This model yields several interesting features not fully explained by 
traditional models: (1) 1/f power spectrum, (2) stable scaling form for the distribution of 
the variations in the beat-to-beat intervals and (3) Fourier phase correlations [f|]. 

To illustrate the approach for the specific example of neuroautonomic control of cardiac 
dynamics, we first note that the healthy heart rate is determined by three major inputs: 

(1) the sinoatrial (SA) node; (ii) the parasympathetic (PS); and (iii) the sympathetic (SS) 
branches of the autonomic nervous system. 

(i) The SA node or pacemaker is responsible for the initiation of each heart beat; in the 
absence of other external stimuli, it is able to maintain a constant interbeat interval [|TJ. 
Experiments in which PS and SS inputs are blocked reveal that the interbeat intervals are 
very regular and average only 0.6s ||. The input from the SA node, Is a, thus biases the 
interbeat interval r toward its intrinsic level t$a (see Fig. |]b). 

(ii) The PS fibers conduct impulses that slow the heart rate. Suppression of SS stimuli, 
while under PS regulation, can result in the increase of the interbeat interval to as much 
as 1.5s ||. The activity of the PS system changes with external stimuli. We model these 
features of the PS input, Ips, by the following conditions: (1) a preferred interval, Tps(n), 
randomly chosen from an uniform distribution with an average value larger than t$a, and 

(2) a correlation time, Tps, during which Tps does not change, where Tps is drawn from a 
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distribution with an average value Ti oc k. 

(hi) The SS fibers conduct impulses that speed up the heart beat. Abolition of parasym- 
pathetic influences when the sympathetic system remains active can decrease the interbeat 
intervals to less than 0.3s ||. There are several centers of sympathetic activity highly sensi- 
tive to environmental influences [jS]. We represent each of the N sympathetic inputs by P ss 
(j = 1, • • • , N). We attribute to I 3 SS the following characteristics: (1) a preferred interbeat 
interval T J ss (n) randomly chosen from a uniform distribution with an average value smaller 
than tsa, and (2) a correlation time Tj in which Tg S (n) does not change; Tj is drawn from 
a distribution with an average value Ti oc k which is the same for all N inputs (and the same 
as for the PS system), so Ti oc k is the characteristic time scale of both the PS and SS inputs. 

The characteristics for the PS and SS inputs correspond to a random walk with stochastic 
feedback control (Fig. |l]c). Thus, for the present example of cardiac neuroautonomic control, 
we have N + 2 inputs and Eq. (0a) becomes 

N 

r(n + 1) - r(n) = I SA (n) + I PS (n, T PS {n)) + £ I 3 3S (n, r J ss (n)) , (4) 

i=i 

where the structure of each input is identical to the one in Eq. (|3|b). Equation (J|) cannot fully 
reflect the complexity of the human cardiac system. However, it provides a general frame- 
work that can easily be extended to include other physiological systems (such as breathing, 
baroreflex control, different locking times for the inputs of the SS and PS systems 0J|, 
etc.). Equation (f|) captures the essential ingredients responsible for a number of important 
statistical properties of the healthy heart rate. 

To qualitatively test the model, we first compare the time series generated by the stochas- 
tic feedback model and the healthy heart and find that both signals display complex vari- 
ability and patchiness (Fig. ^]a). To quantitatively test the model, we compare the statistical 
properties of heart data with the predictions of the model: 

(a) We first test for long-range anti-correlations in the interbeat intervals, which exist for 
healthy heart dynamics ||. These anti-correlations can be uncovered by calculating power 
spectra, and we see (Fig. H) that the model correctly reproduces the observed long-range 
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anti-correlations. In particular, we note that the non-stationarity of both the data and model 
signals leads to the existence of several distinct scaling regimes in the power spectrum of 
r(n) (Figs. and |3[). The stochastic feedback mechanism thus enables us to explain the 
formation of regions (patches) in the time series with different characteristics (see caption 
to Fig. |). 

(b) By studying the power spectrum of the increments we are able to circumvent the 
effects of the non-stationarity. Our results show that l//-scaling is indeed observed for the 
power spectrum of the increments, both for the data and for the model (Fig. |2|). 

(c) We calculate the probability density P(A) of the amplitudes A of the variations of 
interbeat intervals through the wavelet transform. It has been shown that the analysis of 
sequences of interbeat intervals with the wavelet transform can reveal important scaling 
properties || for the distributions of the variations in complex nonstationary signals. In 
agreement with the results of Ref. ||, we find that the distribution P(A) of the amplitudes 
A of interbeat interval variations for the model decays exponentially — as is observed for 
healthy heart dynamics (Fig. ^). We hypothesize that this decay arises from nonlinear 
Fourier phase interactions and is related to the underlying nonlinear dynamics. To test this 
hypothesis, we perform a parallel analysis on a surrogate time series obtained by preserving 
the power spectrum but randomizing the Fourier phases of a signal generated by the model 
(Fig. P(A) now follows the Rayleigh distribution P(A) ~ Ae~ A2 [10), since there are no 
Fourier phase correlations. 

(d) For the distribution displayed in Fig. |J we test the stability of the scaling form at 
different time scales; we find that P(A) for the model displays a scaling form stable over a 
range of time scales identical to the range for the data (Fig. |) || . Such time scale invariance 



indicates statistical self-similarity [11 



A notable feature of the present model is that in addition to the power spectra, it accounts 
for the functional form and scaling properties of P(A), which as we show are independent 
from the power spectra. 

The model has a number of parameters, whose values may vary from one individual to 



another, so we next study the sensitivity of our results to variations in these parameters. We 
find that the model is robust to parameter changes. The value of T lock and the strength of 
the noise a are crucial to generate dynamics with scaling properties similar to those found for 
empirical data. We find that the model reproduces key features of the healthy heart dynamics 
for a wide range of time scales (500 < T lock < 2000) and noise strengths (0.4 < a < 0.6). 
The model is consistent with the existence of an extended "zone" in parameter space where 
scaling behavior holds, and our picture is supported by the variability in the parameters for 
healthy individuals for which similar scaling properties are observed. 

The model, and the data which it fits, support a revised view of homeostasis that takes 
into account the fact that healthy systems under basal conditions, while being continuously 
driven away from extreme values, do not settle down to a constant output. Rather, a more 
realistic picture may involve nonlinear stochastic feedback mechanisms driving the system. 



6 



REFERENCES 

[1] C. Bernard, Les Phenomenes de la Vie (Paris, 1878); W. B. Cannon, Physiol. Rev. 9, 
399 (1929). 

[2] M. F. Shlesinger and B. J. West, in Random Fluctuations and Pattern Growth (Kluwer, 
Dordrecht, 1988). 

[3] G. H. Weiss, Aspects and Applications of the Random Walk (Elsevier Science B.V., 
North- Holland, New- York, 1994). 

[4] M. G. Rosenblum and J. Kurths, Physica A 215, 439 (1995); H. Seidel and H. P. Herzel, 
in Modelling the Dynamics of Biological Systems E. Mosekilde and O. G. Mouritsen, 
Eds. ( Springer- Verlag, Berlin, 1995); J. M. Hausdorff and C. -K. Peng, Phys. Rev. E 
54, 2154 (1996). 

[5] R. M. Berne and M. N. Levy, Cardiovascular Physiology 6th ed. (C.V. Mosby, St. Louis, 
1996). 

[6] C. -K. Peng et al, Phys. Rev. Lett. 70, 1343 (1993). 

[7] I. Daubechies, Comm. Pure and Appl. Math. 41, 909 (1988). 

[8] J. F. Muzy, E. Bacry, and A. Arneodo, Int. J. Bifurc. Chaos 4, 245 (1994); A. Arneodo et 
al, Physica D 96, 291 (1996). 

[9] P. Ch. Ivanov et al, Nature 383, 323 (1996). 

[10] R. L. Stratonovich, Topics in the Theory of Random Noise (Gordon and Breach, New- 
York, 1981). 

[11] J. B. Bassingthwaighte, L. S. Liebovitch, and B. J. West, Fractal Physiology (Oxford 
Univ. Press, New York, 1994). 



7 



FIGURES 




8 



FIG. 1. Schematic representation of the dynamics of the model, (a) Evolution of a random 
walk starting from initial position ro- The deviation of the walk from level ro increases as n 1 / 2 , 
where n is the number of steps. The power spectrum of the random walk scales as l// 2 (Brownian 
noise). The distribution P(A) of the amplitudes A of the variations in the interbeat intervals 
follows a Rayleigh distribution, (b) Random walk with a bias toward ro- For short time scales 
(high frequencies), the power spectrum scales as l/f 2 (Brownian noise) with a crossover to white 
noise at longer time scales due to the attraction to level ro . Note the shift of the crossover to longer 
time scales (lower frequencies) when stronger noise is present. However, in both cases, P(A) follows 
the Rayleigh distribution, (c) Random walk with two stochastic feedback controls. In contrast to 
(b), the levels of attraction ro and t\ change values in time. Each level persists for a time interval 
Tj drawn from a distribution with an average value Ti oc k- Perturbed by changing external stimuli, 
the system nevertheless remains within the bounds defined by At even after many steps. We find 
that such dynamical mechanism, based on a single characteristic time scale ?i oc k, generates a l/f 
power spectrum over several decades. Moreover, P(A) decays exponentially, which we attribute to 
nonlinear Fourier phase interactions in the walk. 
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FIG. 2. (a) Sequence of interbeat intervals r for a healthy individual, (b) Sequence of 
interbeat intervals for the model with parameters N = 7, wsa = wss = wps/3 = O.Olsec. For 
the results presented here, Tj is randomly chosen from an exponential distribution with average 
T\ock = 1000 beats, and rj is a symmetrical exponential distribution with zero average and standard 
deviation a = 0.5. The preferred values of the interbeat intervals for the different inputs were 
picked according to the following rules: (1) tsa = 0.6sec, (2) tps are randomly selected from an 
uniform distribution in the interval [0.9, 1.5]sec, and (3) the r^ 5 's are randomly selected from an 
uniform distribution in the interval [0.2, 1.0]sec. We note that the actual value of the preferred 
interbeat intervals of the different inputs and the ratio between their weights are physiologically 
justified and are of no significance for the dynamics — they just set the range for the fluctuations of 
r, chosen to correspond to the empirical data. Also, any change of the shape of the distribution of 
the noise term or of the locking times leaves the reported here statistical properties of the generated 
signal unchanged, (c) Power spectra of the interbeat intervals r(n) from the data and the model 
described by the relation S(f) ~ l// 1 " 1 . The presence of patches in both heart and model signals 
lead to observable crossovers embedded on this 1/f behavior at different time scales. The local 
exponent (3 from the power spectrum of 24h records (« 10 5 beats) for 20 healthy subjects shows a 
persistent drift, so no true scaling exists, (d) Power spectra of the increments in r(n). The model 
and the data both scale as power laws with exponents close to one. Since the non-stationarity is 
reduced, crossovers are no longer present. Here the local exponent (3i fluctuates around an average 
value close to one, so true scaling does exist. 
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FIG. 3. Top: Effect of the correlation time Ti oc k on the scaling of the power spectrum of r(n) 
for a generated signal comprising 10 6 beats. With increasing Ti oc k, the power spectrum does not 
follow a single power law but actually crosses over from a behavior of the type l/f 2 at very small 
time scales (or high frequencies), to a behavior of the type 1//° for intermediate time scales, 
followed by a new regime with l/f 2 for larger time scales. At very large time scales, another 
regime appears with flat power spectrum. Bottom: Schematic diagram illustrating the origin of 
the different scaling regimes in the power spectrum of r(n). For very short time scales, the noise will 
dominate, leading to a simple random walk behavior and l/f 2 scaling (Region A). For time scales 
longer than T4, the attraction towards the "average preferred level" of all inputs will dominate, 
leading to a flat power spectrum (Region B, see also Fig. |l|b). However, after a time Tb (of the 
order of T\ oc ^/N), the preferred level of one of the inputs will have changed, leading to the random 
drift of the average preferred level and the consequent drift of the walker towards it. So, at these 
time scales, the system can again be described as a simple random walker and we expect a power 
spectrum of the type l/f 2 (Region C). Finally, for time scales larger than Tc, the walker will start 
to feel the presence of the bounds on the fluctuations of the preferred levels of the inputs. Thus, 
the power spectrum will again become flat (Region D). 
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FIG. 4. We apply to the signal generated by the model the wavelet transform with fixed scale 
a, then use the Hilbert transform to calculate the amplitude A. Top: normalized histogram P(A) 
for the data (6h daytime) and for the model (with the same parameter values as in Fig. ||), and 
for a = 8 beats (~ 40s). The inset shows a similar plot for data collected during the night and 
for the model with N < wps/wss- Note that the distribution is broader for this case with large 
values for the amplitudes deviating from the exponential tail. Bottom: we test the stability of 
the analysis for the model at different time scales a. The distribution is stable over a wide range 
of scales (identical to the range observed for heart data) and indicates statistical self-similarity in 
the variations at different time scales. We test the generated signal for nonlinearity and Fourier 
phase correlations, creating a surrogate signal by randomizing the Fourier phases of the generated 
signal but preserving the power spectrum (thus, leaving the results of Fig. [2] unchanged). The 
histogram of the amplitudes of variations for the surrogate signal follows the Rayleigh distribution, 
as expected theoretically (see inset). Thus the observed distribution which is universal for healthy 
cardiac dynamics, and reproduced by the model, reflects the Fourier phase interactions. 
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